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Abstract 

In this work, we develop and apply the multi-layer multi-configuration time-dependent Hartree 
method for bosons (ML-MCTDHB), which represents a unique tool for investigating the non- 
equilibrium dynamics of multi-species bosonic systems in arbitrary dimensions. Being an ab initio 
method for solving the time-dependent Schrodinger equation, ML-MCTDHB takes all correlations 
into account. The multi-layer feature of ML-MCTDHB allows for tailoring the wave function 
ansatz in order to describe intra- and inter-species correlations accurately and efficiently To show 
the beneficial scaling and the efficiency of the method, we explore the correlated dynamics of three 
species tunneling in a double well trap. Here we demonstrate and analyze in detail the build up 
of inter- and intra-species correlations in the course of the quantum dynamics as well as signatures 
of equilibration. Furthermore, we show that the presence of the third species which ultra-weakly 
attracts bosons of the two other species can reduce the correlations between those two species. 
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I. INTRODUCTION 



Due to the high degree of controllability and isolatedness, ensembles of ultracold atoms 
have been attracting much attention for almost two decades. Whereas the first experiments 
with ultracold bosons were situated in parameter regimes which allow for a mean field de- 

nn 

scription by means of the famous Gross Pitaevskii equation [lH3|, nowadays experimental 
techniques enable us to also study the effects of strong interactions leading to highly corre- 
lated states \4\. Furthermore, different trapping geometries resulting in lower dimensional 
ensembles as well as effects due to the presence of other components (e.g. are in the 
focus. In particular, it has become feasible to trap and manipulate ensembles of only a few 
atoms (e.g. s|, which allows to explore the impact of the particle number on various 
effects. The quest for a description of the non-equilibrium quantum dynamics beyond the 
mean field approach, i.e. by taking many-body correlations into account, demands new 
theoretical tools. 

Simulating the quantum dynamics of an interacting iV particle system, however, is a 
tough task in general due to the exponential scaling of the state space: If, for instance, 
one allows each particle to occupy M single particle states, the many-body Hilbert space 
turns out to be M N dimensional for a system of distinguishable particles. When going to 
higher particle numbers, one, hence, has to truncate the Hilbert space further. Moreover, 
this scaling problem has a serious consequence for simulating dynamics: In the course of 
time, the system state will move through different subspaces of the total Hilbert space in 
general. If one now tries to expand the system state with respect to a time-independent 
basis, many basis states will be needed in order to span the relevant subspace for the whole 
propagation time. So for reducing the number of basis vectors, one should better employ 
a time-dependent, moving basis, in which the instantaneous system state can be optimally 
represented. 

In the multi-configuration time-dependent Hartree method (MCTDH), such a co-moving 
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basis is provided by time-dependent Hartree products [8, 9]. The equations of motion for 
the single particle functions (SPFs) building up the time- dependent Hartree products as 
well as for the expansion coefficients are obtained by means of the Dirac-Frenkel variation 
principle, which ensures a variationally optimal representation of the total wave function at 
any instant. Using Hartree products as the many-body basis states, MCTDH is well suited 



for dealing with distinguishable particles. Originally developed in the context of quantum 
molecular dynamics, MCTDH has successfully been applied to an enormous diversity of 
problems concerning the vibrational and wave packet dynamics of molecules (cf. Jiol . ll| 
and reference therein). Moreover, MCTDH turned out to be also powerful for simulating 
strongly correlated few-body systems of ultracold, indistinguishable bosonic atoms in traps, 
see e.g. [12 



MCTDH has been extended for dealing with higher particle numbers in several ways: In 
multi-layer MCTDH (ML-MCTDH), degrees of freedom are combined to higher-dimensional 
ones. The MCTDH expansion scheme is then applied to these combined modes and their 
constituting degrees of freedom in a cascade. The ML-MCTDH theory has been presented 
and implemented for the first time in [141 ] . Afterwards, a recursive formulation of the ML- 
MCTDH theory was derived, which allows for a fast and convenient treatment of arbitrarily 
complex ML-MCTDH expansion schemes 15|. The ML-MCTDH method is particularly 



fruitful for scenarios where there are strong correlations within certain subsystems while 
the inter-subsystem correlations are relatively weak 1614221] . Especially, many system-bath 
problems are of this type. In such a situation, the necessary number of time-dependent 
basis functions can be reduced by combining the strongly correlated degrees of freedom in 
an ML-MCTDH expansion. As MCTDH, ML-MCTDH is based on expansions in terms of 
Hartree products and, thus, ignores possible particle exchange symmetries of the system. 

In contrast to this, the MCDTH method for fermions (MCTDHF) and for bosons (MCT- 
DHB) reduce the number of expansion coefficients by expanding the total wave function 



in terms of slater determinants and permanents 
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-26], respectively. MCTDHF has been 



applied in the context of quantum chemistry as well as for molecules in strong laser fields 
27H30]. The crossover from the Gross-Pitaevskii to the many-body regime has been studied 



for several ultracold bosonic systems by means of MCTDHB (e.g. 3l|, |32|) . There is even 



a formalism, which unifies MCTDHB and MCTDHF for dealing with bose-bose, bose-fermi 



and fermi-fermi mixtures 33] . Furthermore, an extension for considering particle conversion 
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in a multi-species system has been developed [34]. 

ML-MCTDH features a high flexibility with respect to the wave function expansion: Any 
way of combining degrees of freedom to higher dimensional ones in a cascade is allowed and 
leads to a specific wave function ansatz. As, in the case of indistinguishable particles, many 
of these mode-combination schemes do not reflect the indistinguishability of the particles, 



the particle exchange (anti-) symmetry can, in general, not be incorporated in ML-MCTDH, 
making ML-MCTDH incompatible with MCTDHB and MCTDHF. There is, however, an 
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exception: In 
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DH- 



37|, for 



35] . a second quantization representation of ML-MCDTH (M. 
SQR) has been derived, which has been applied to charge transport problems 
example. The ML-MCTDH-SQR method is based on the factorization of the Fock space 
for m single particle modes into m Fock spaces corresponding to just one single particle 
mode each, which leads to the product structure required by ML-MCTDH. In this way, 
ML-MCTDH-SQR benefits from the flexible cascade scheme of ML-MCTDH while keeping 
the SPF, constituting the Fock states, time-independent. 

For investigating complex bosonic systems, a new ab initio method, called multi-layer 
multi-configuration time-dependent Hartree method for bosons (ML-MCTDHB), is derived 
and applied to explore novel quantum dynamics in this work. This method provides a new 
symbiosis of MCTDHB and ML-MCTDH, which in particular differs from ML-MCTDH-SQR 
in the way that it is based on time-dependent SPFs. To face the challenge of simulating the 
non-equilibrium dynamics of complex bosonic systems without discarding any correlations, 
we derive two basic classes of ML-MCTDHB wave function ansatzes: (i) a multi-layer ex- 
pansion with respect to the different bosonic species for optimally taking account inter- and 
intra-species correlations and (ii) a particle multi-layer expansion in which the single particle 
functions are further expanded in the MCTDH form. The latter allows for an optimal wave 
function ansatz when dealing with trap geometries in higher dimensions or internal degrees 
of freedom. In both classes of expansion schemes, the bosonic exchange symmetry is exactly 
taken into account, which improves the accuracy as well as the efficiency of the method. 
Furthermore, any combination of the expansion schemes (i) and (ii) is possible resulting in 
a highly flexible ab initio approach to systems with - in principle - arbitrary many species 
in arbitrary dimensions. 

This paper is organized as follows: Firstly, we will review the methods ML-MCTDH 
and MCTDHB in sections [HA] and [HB] The wave function expansion of ML-MCTDHB 
is consecutively introduced in section III C II and afterwards the derivation principle for the 
equations of motion as well as the results for multi- and single-species systems are presented 
in sections CTT721 lH"C3l and [HUH respectively. In section [HQ ML-MCTDHB is applied 
to several tunneling scenarios, which show interesting correlation effects and at the same 
time illustrate how the ML-MCTDHB wave function ansatz can be adapted to capture 



the relevant physics. We conclude with a summary, discussing also an extension of ML- 
MCTDHB to fermionic species in section HVl Finally, we would like to note that the following 
reviews of the existing methods ML-MCTDH and MCTDHB are indispensable in order to 
be self-contained, to embed ML-MCTDHB into the existing methods, to give a consistent 
definition of the layer structure in connection with particle statistics and in order to be able 
to present ML-MCTDHB in its full generality. 



II. THEORY 



A. ML-MCTDH 



Firstly, let us summarize the MCTDH and ML-MCTDH approaches [SHlH, H4Hl6|, |38||. 
We will review these two methods in such a way that an extension to the subject of this work, 
ML-MCTDHB, can be achieved smoothly. At the first place, we will introduce the wave 
function expansion schemes for both MCTDH and ML-MCTDH. Afterwards, the equations 
of motion are discussed. 

Let us consider the general situation of a system consisting of N physical degrees of 
freedom (DOF) (xi, x 2 , x^), where the DOF x K can correspond to a position, momentum 
or internal degree of freedom. Furthermore, N does not need to correspond to the number 
of physical particles but can take internal degrees of freedom as well as the dimensionality of 
the system into account. For instance, a system of N physical particles in three-dimensional 
space can either be modeled by 3N DOF, (xi,X2, ...,Xsn), or by N three-dimensional DOF, 
(xi, X2, ...,xn)- To each physical DOF k, we then associate a set of orthonormal time- 
dependent single particle functions (SPFs) {<f)^(x K ),n = l,...,m K }. In MCTDH, the state 
vector of the system is expanded with respect to the (truncated) basis of Hartree products 
constituted by the SPFs of each DOF: 

i*(*)>= E ^«(*)i^(*))i^(*))---iC(*))^E^i J )- w 
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To step from MCTDH to ML-MCTDH, we first perform a mode combination of the N 
physical DOF and group them into Np_i so-called logical DOF such that the first n\ physical 
DOF are combined to the logical DOF a;f _1 , the next n 2 physical DOF are combined to the 



logical DOF x% 1 , and so on: 
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X-y (x\, X2, x ni ), 
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P-1 _ P-1 ( \ 
X N P _ 1 ~ X N P _ 1 \ x N-n Np + l, XN-n Np +2, Xn)- 

The physical DOF {xi} are said to form the physical layer in ML-MCTDH, and the logical 
DOF {x { _1 } constitute the first upper layer above the physical layer, where the superscript 
P — 1 denotes that these DOF are one layer above the physical layer. We can continue 
and combine now the iVp_i logical DOF {xf -1 } to generate a new set of Np_ 2 (Np^ 2 < 
Np-i) logical DOF {xf ~ 2 | % G [l,Np- 2 ]}, which form the second upper layer above the 
physical layer. By continuing this mode-combination scheme until reaching a layer containing 
only one logical DOF, we complete the construction of the so-called tree structure in ML- 
MCTDH. Generally speaking, the tree structure of ML-MCTDH contains different layers of 
DOF, and the DOF on each layer a are generated via mode combination of the corresponding 
DOF on the lower lying layer a + 1. 

Similarly to the physical DOF, we also associate a set of orthonormal time-dependent 
SPFs {\<Pn' R )} to each logical DOF a£. Please note that, despite of the term "SPF", the state 
\(f>n' K ) can correspond e.g. only to the x- and y-coordinate of a particle in three-dimensional 
coordinate space or to the x-coordinates of four different particles, depending on the mode- 
combination scheme. 

Each of these SPFs is expanded with respect to Hartree products made of the SPFs of 
those DOF which constitute the logical DOF x a - 

10= E <t«,..,i m iC l! " 1 >i^ lw >---iC w > = E^i itt;,i >- ( 3 ) 
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We note that the time-dependencies are omitted in the notation. Here and from now on, 
the ingredients of all formulas have to be regarded as being time-dependent, unless stated 
otherwise. 

Especially, the single DOF of the top layer possesses only a single SPF, which turns out 
to be the state vector |\&) of the system and is expanded in terms of Hartree products made 
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of the SPFs of the second layer: 

i*> = E^i j1 >- ( 4 ) 
i 

Equations ([3]) and (TJj) indicate that each SPF is "locally" expanded in the MCTDH manner 
(PQ), in terms of Hartree products made of the SPFs of the lower layer. Finally, the ML- 
MCTDH wave function representation is completed by expanding the SPFs of the physical 
DOF with respect to some given time-independent basis function \x K ), which one usually 
calls primitive basis functions (PBFs): 



J-K 



The physical DOF x K can be high- dimensional, for instance, x K might contain the x-, y-, and 
z-position degrees of freedom for the case of particles in three dimensional coordinate space. 
In this example, the SPFs of the DOF x K would be expanded in terms of three-dimensional 
PBFs \x K ) = \r K ), which is known as mode combination in MCTDH (cf. |8|). We will call 
such a mode combination primitive mode combination in order to distinguish it from the 
mode combinations in the other layers of a ML-MCTDH wave function expansion. Finally 
we remark that having the sets of PBFs fixed, the evolution of the system turns out to 
be completely represented by the time evolution of the expansion coefficients of each DOF, 

i\I 

Before we state the equations of motions for those coefficients A^f, we briefly review 
the diagrammatic representation of the wave function expansion ( l3]Hf5j) introduced by 
Manthe 15]: Performing the MCTDH expansion layer by layer in a cascade naturally leads 
to a visualization in terms of a tree diagram. In such a tree diagram, every node is identified 
with a logical or physical DOF x™. Its number of teeth represents how many SPFs \(f\ ,K ) 
are provided for this DOF. Then each link between two nodes indicates the summation 
over an index in the expansion (l3ll4l[5|) . Finally, a rectangular box stands for a set of PBFs 
and the number of its teeth coincides with the number of given PBFs. Figure [1] illustrates 
this diagrammatic notation by an example of two particles in two-dimensional coordinate 
space: The physical layer consists of four DOF corresponding to the two spatial directions 
for each particle. In the second layer, two DOF each are combined to one logical DOF each, 
representing the two physical particles. In the following we will equivalently refer to the 
DOF as nodes. This will enable us to address the DOF, which turn out to be relevant for 



a node under consideration, by the pictorial terms parent and child node, respectively. At 
this step, it is clear that MCTDH could be seen as a special case of ML-MCTDH, which 
contains only two layers: the top layer of and the primitive layer of physical DOF. 

In this way, we will mainly focus on the equations of motion of ML-MCTDH, and those of 
MCTDH are naturally included. 




Figure 1: Diagrammatic representation of a ML-MCTDH wave function expansion for a 

system of two physical particles, which correspond to the nodes of the second layer. 
Indicated by their teeth, four and two SPFs are provided for the first and second physical 

particle, respectively. According to the third layer, each of these SPFs is unraveled by 
Hartree products of two one-dimensional SPFs, i.e. orbitals of the x- and y-direction. Each 
of these one-dimensional orbitals is in turn represented by n x and n y , respectively, PBFs 
or, pictorially speaking, spatial gridpoints. For illustrating the notation, the expansion 
coefficients are written into the corresponding nodes. 

We assume the Hamiltonian of the system to be given as a sum over products of single par- 
ticle operators acting on each physical DOF, which is particularly suitable for ML-MCTDH 
as it allows for a recursive formulation of the theory layer by layer: 

JV 

V K=l V 

Here, denotes a single particle operatoJl] acting on the physical DOF x K . Some interaction 
1 Please note that the term "single particle operator" does not need to refer to a physical particle but refers 
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potentials can be brought analytically into the desired product form see the harmonic 
interaction v(x\ — 22) oc (x\ — X2) 2 , for example. Moreover, such a decomposition of a 

ized by 



39 



40J. 



many-body operator into products of single particle operators can always be rea 
the so-called POTFIT algorithm, which essentially solves iV eigenvalue problems 
For a Hamiltonian with up to F— body operators, the number of v terms in in general 
equals the product of the dimensions of the Hilbert spaces corresponding to the respective 
F physical DOF but one of these. The v term summation, however, can be truncated in a 
controlled way, and often, only a few v terms are necessary for a fair representation of the 
total Hamiltonian. In this case, the product form ([!)]) can improve the efficiency because 
high-dimensional integrals decompose into low-dimensional ones when calculating matrix 
elements. We remark that the Hamiltonian H may be time-dependent in ML-MCDTH. 

Since every layer in itself provides a complete description of the system, the Hamiltonian 
can be also described equivalently by the SPFs of each layer: For the a layer, the single 
particle operators can be grouped according to the mode-combination scheme up to the 
layer a, which leads to the a— layer representation of the Hamiltonian H: 

H = J2 c »U h » K = J2 c ^:- (7) 

V K=l V 

Here, h"' ,K denotes an operator exclusively acting on the logical DOF x®. Hence, it can be 
expressed in terms of the SPFs of this DOF as well as the operators of the layer a + 1, whose 
product equals h"' Kl : 

K K )nm = = E AZil^l II K +li V aiK ), (8) 

I, J i=K\ 

where the DOF (a; k) is a combination of the DOF (a + 1, Ki) to (a + 1, n n ). 

A useful concept for the equations of motion is the single hole function (SHF) of the node 
(a; k), which is given by the projection of the total wave function onto the i-th SPF 

of the node (a; k): 

i*r> = (c ;K i*>- (9) 

The SHF of a node (a; k) can be calculated by expanding the total wave function in terms of 
the SPFs of the layer a and performing the projection ([9]). A recursive formula for obtaining 



the SHFs layer by layer is stated in 
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to a physical DOF. 
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Now we can substitute the ansatz of \^) to the Dirac-Frenkel variational principle 4ll.l42|: 



(6*\(iat-H)\*) = 0, (10) 
to obtain the equations of motion for the coefficients of each layer. The Dirac-Frenkel 



variationa 
one 
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principle turns out to be equivalent to the Lagrangian and to the McLachlan's 
44j. The latter essentially says: Propagate your ansatz \^f) according to idt\^f) = |@) 
with |©) being of the ansatz form and minimizing \\i\Q) — H\^)\\ 2 . This ensures that we 
obtain a variationally optimal wave function within our class of ansatzes. Let us turn back 
to (ITUj) . which is easier to handle: The variation with respect to the top node coefficients, i.e. 
6(1))* (^\(idt — H)\i&) = 0, effectively gives the projection of the Schrodinger equation on the 
Hartree product {I 1 ). This results in the equation of motion for the top layer coefficients: 

l d t A] = Y.^\H\J 1 )A 1 J . (11) 
j 

The variation with respect to a non-top node coefficient results in effectively projecting 
the Schrodinger equation on the SHF times the Hartree product |/ Q;K ). After some 

algebra, one obtains the equation of motion for the coefficients A"'*: 

id t A%« = ^a{^(/ a; 1(l - P a ' K )hT\J a '' K )}* (12) 

V J 

m _ 

Here, \{p a ' K )~ l * (H^' ,K )] denotes the matrix product of the inverscl of the density matrix of 
DOF i° and the so-called mean field matrix of the v-th Hamiltonian term with respect to 
x": The density matrix p a ' ,K is defined as the overlap of different SHFs of the DOF x°: 

It turns out that equals the transposed of the physical reduced density matrix of the 
DOF x®, which one obtains by tracing out all DOF but x" in the many-body state 
The mean field matrix (if" ;ft ) n m is defined as the scalar factor of the mean field operator 
when making use of the product structure of (cf. (|7j)): 



As a density matrix can be singular and will be singular especially for converged MCTDH type simulations, 



one has to regularize the density matrix in (|12p . which can be done without affecting the total wave function 

(cf. 3). 
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which can be evaluated by means of OH]) in a recursive manner [15]. The operator P a]K equals 
the projection operator onto the space spanned by the SPFs of the node (a; k). 

We note the following interpretation of the terms in (fT2[) : The operator h"' ,K is responsible 
for the transition between different Hartree products \J a ' ,K ), whereas the density and the 
mean field matrices scatter the SPF to |0" ;ft ). Furthermore, it has to be pointed out 

that the equations of motion f lTTfT2|) are valid under the additional set of constraints: 

(d t <f>T K W K ) = °> ( 15 ) 

which obviously ensures the orthonormality of the SPFs. In the following sections, we 
implicitly impose the constraint ([TBI) on all occurring SPFs since such an additional condition 



18 



is needed for obtaining unique equations of motions (cf. 

Finally, we would like to emphasize that the flexible cascade of MCTDH expansions 
([3]fT|B]) leads to an optimal resolution of the relevant subspace of the many-body Hilbert 
space within the truncated basis given by the SPFs. Nevertheless, both MCTDH and ML- 
MCTDH are methods tailored for treating distinguishable DOF. In the situation of identical 
particles, i.e. fermions or bosons, the Hilbert space is restricted by the permutation (anti-) 
symmetry. Applying MCTDH directlju to such a system is rather inefficient due to the 
redundancies in the expansion coefficients. However, it is possible to modify MCTDH by 
taking into account the permutation symmetries, which can significantly improve the effi- 
ciency In the following section we will focus on the treatment of bosonic systems and briefly 
introduce the MCTDHB, which specializes MCTDH to the bosonic system for achieving 
higher performance. 

B. MCTDHB 

We now consider a single-species system of N bosons. All bosons share the same set 



23|-[25|, the many- 



of time-dependent orthonormal orbitals {|0i(£))}ie[i,Afl- 111 MCTDHB 
body wave function is then expanded within the basis of permanents, i.e. bosonic number 
states \n;t), which naturally takes the bosonic symmetry into account and, thus, reduces 
the number of expansion coefficients in comparison to a MCDTH expansion. Omitting all 



3 Since the Hamiltonian for a system of identical particles preserves the (anti-) symmetry of the total wave 
function during propagation, so does a converged MCTDH propagation. 
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time- dependencies in the notation as always, this yields the following ansatz wave function: 

\^} = J2Cn\n} } (16) 

n\N 

where n = (m, ...,um) and rti is the occupation number of bosons in the orbital \<fii). The 
symbol "n\N" shall indicate that the summation runs over all positive integers with the 
constraint Y2fLi n i = N. The number states |n) become time-dependent as they equal 
permanents of time-dependent orbitals: 

|n> = (AT! ni !...n M !)-i ^ ' -\K {n) )n- (17) 

This equation is based on a multi-index I = (zi, ...,2at) which is fixed such that j occurs rij 
times for any j G [1,M]. The summation then runs over all permutations ir of the first N 
integers, i.e. ir : [1,N] — > [1,N], and ket vectors get another index, i.e. | ) K , encoding to 
"which" boson the state belongs to. Furthermore, the time-dependent orbitals are in turn 
expanded with respect to some time- independent PBFs {|x)}: 

i<^> = ]r<M*>. (is) 

X 

Next, we extend the diagrammatic representation of the total wave function to expansions 
in terms of permanents, which will turn out to be useful when discussing the expansion 
schemes of ML-MCTDHB: In the following, we indicate with a plus in a node that its SPFg^ 
are expanded in permanents. Then, the number of vertical lines below that bosonic node 
represents the number of bosonic particles, while the number of teeth of its child node equals 
the number of single particle orbitals. An example for a diagrammatic representation of a 
MCTDHB wave function for a system of six bosons in one- dimensional coordinate space is 
given in figure [2j 

When restricting to at most two-body interactions, the many-body Hamiltonian can be 
expressed in the second-quantization picture as follows: 

H = h k q a{a q + - ^2 v ksi q a\a\a l a q (19) 

k,q k,s,l,q 

where a\{a,k) is the time-dependent bosonic creation (annihilation) operator corresponding to 
the time-dependent orbital \(f)f.) and obeying [a s ,a g ] = 0, [a s ,aj] = 5 sq . The time- dependent 



4 Note that in a MCTDHB tree, the top node is always bosonic and represents a single SPF, namely the 
total wave function. 
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N = 6 



m = 3 



71 = 



Figure 2: Visualization of a MCTDHB wave function expansion for a system of iV = 6 
bosons, which may occupy m = 3 orbitals. Each orbital is represented by n = 5 

time-independent PBFs. 



matrix elements of the one-body Hamiltonian h and two-body interaction potential v, i.e. 
h^q and Vksiq respectively, are defined as: 

h kq = ((f) k \h\(f) q }, 

Vkslq = (4>k\(4>s\v\4>l}\4>q}, 

and have to be evaluated in the representation f fT8l) . 

The equations of motion for and {|0fc(^))} can be derived from the Dirac-Frenkel 
variational principle and read: 



id t C H = y^(n|g|m)C„- 



rh\N 



M 



id, 



t\<Pk 



> = (1 - P) h\<f>k) + ^2 [P %j PjslqVsl\<Pq) 



(21) 
(22) 



In the latter equation, P is the projection operator onto the subspace spanned by the 
SPFs Pjk denotes the element of the so-called one-body density matrix, defined as 

Pjk = (ty\a^cik\ty) , while pksiq is the element of the so-called two-body density matrix, defined 
a s Pkslq = (^\aldldid q \^). Finally, the one-body Hamiltonian h as well as the mean field 
operator v s i can be expressed in terms of the PBFs as: 



v s i 



h= } j (x 1 \h\x 2 ) |xi)(x 2 |, 

Xl,X2 

22{xi\{(f) s \v\(f)l)\x2) l^lX^I- 



(23) 
(24) 



X1,X2 
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By projecting (122]) on a PBF \x), one arrives at the equation of motion for the coefficients 
4>i(x). We note that (122]) is equivalent to the ML-MCTDH non-top node equations of motion 
([12]) . The special form of (122]) is a consequence of the Hamiltonian ( [19]) containing only up to 
two body operators. The approach sketched above can be generalized to bose-bose mixtures 



and bose-fermi mixtures, and we refer the reader to 



24 



25] for further details. 



C. ML-MCTDHB 

ML-MCTDH has proven to be a powerful tool for simulating complex systems with many 
degrees of freedom, while MCTDHB is optimized for treating scenarios involving one bosonic 
species. In this chapter, we derive the novel ML-MCTDHB method which both explicitly 
considers the bosonic symmetry of the particles and allows for a flexible adaption to system 
specific correlations and geometric aspects. 



1. Concepts of ML-MCTDHB 

While the wave function decomposition of ML-MCTDH is arbitrary in principle, ML- 
MCTDHB is based on two classes of wave function expansions: (i) species multi-layer and 
(ii) particle multi-layer MCTDHB. A general ML-MCTDHB wave function expansion then 
turns out to be a combination of these two classes. In the following, we discuss these classes 
before deriving the equations of motions: 

(i) Species multi-layer MCTDHB: If the system under consideration consists of 5* bosonic 
species, S > 1, the total wave function is expanded in Hartree products Yl K =i l^i'*)' 

Mi M s S 

i*>=e---e<..^iiic>- 

ii = 1 ig=l k=1 

In the following, we refer to \4>f K ) as the i-th species SPF of species k. Hence, the expansion 
in terms of Hartree products reflects the distinguishability of the species. Each species SPF 
\(j)f K ) is in turn represented as a superposition of permanents |n), paying respect to the 
bosonic symmetry: 

I0- ;K > = E4^>- ( 26 ) 

n\N K 
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Here, n = (nj, n mK ) denotes the occupation numbers of the m K single particle orbitals, 
which can be populated by the bosons of species k. The notation n\N K indicates the con- 
straint that the occupation numbers rij have to sum up to the number of bosons of species 
k, i.e. N K . Finally, the SPFs underlying the time-dependent permanents |n) are denoted 
by \($' K ), i G {1, ...,m K }, and expanded with respect to time-independent basis functions as 
usual for the SPFs of a primitive layer in ML-MCTDH. In the following, this class of wave 
function decomposition is denoted as species multi-layer. 

The strength of species multi-layer MCTDHB lies in its flexibility with respect to the 
species SPFs: If the intra-species correlations are much stronger than the inter-species ones, 
only a few species many-body states, which are variationally optimally rotated by the ML- 
MCTDHB equations of motion, might be sufficient for a fair representation of the total 
state. So depending on the details of the system, species multi-layer MCTDHB allows to 
choose the optimal number of species SPFs, ranging from M — 1, the mean field limit 

for weakly interacting species, up to the full Cl-limit for strongly correlated species when 

M= (N +m -iy 
V m— 1 / 

A general species multi-layer expansion for a system of S species, each of them consisting 
of N K bosons, is shown in figure [3j For each species k G {1, if>}, M K species SPFs are 
provided and the N K bosons constituting species k may occupy m K orbitals. Although 
theses m K orbitals are expanded in n K time-independent states in figure [31 it should be 
pointed out that primitive mode combination can also be applied on this layer for a natural 
treatment of two- and three-dimensional bosons or internal degrees of freedom. 

(ii) Particle multi-layer MCTDHB: A multi-layer expansion can also be applied to the or- 
bitals of the bosons if these are higher-dimensional objects living in two- or three-dimensional 
coordinate space and/or having internal degrees of freedom. This multi-layer extension al- 
lows for grouping strongly correlated particle degrees of freedom together leading to a tailored 
representation of the bosonic orbitals for a given system. 

Figure H] shows some examples for particle multi-layer MCTDHB: The dashed lines head- 
ing to the bosonic nodes shall indicate that the shown tree diagrams may be incorporated 
in a species multi-layer MCTDHB expansion. Particularly, we would like to highlight how 
quasi-one or -two dimensional systems can be treated in quite a natural way: Whereas all 



three spatial degrees of freedom are treated equivalently in figure 4(a) , there is an inter- 



mediate layer consisting of the combined xy and the z degrees of freedom in figures 4(c" 
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Figure 3: Archetype of a species multi-layer MCTDHB expansion of the wave function 
expansion for a system of S species, each of them consisting of N K bosons, which may 
occupy m K orbitals, which in turn are expanded in n K time- independent basis function. 
Each species may populate M K species SPFs. Please note that the vertical and horizontal 
lines are only schematic because the numbers of bosons, SPFs and PBFs are arbitrary. As 
S is arbitrary, the dashed lines shall indicate the expansion of the species k' with 

1 < k' < k or k < k' < S. 



and 4(d) This allows for associating much fewer states, here only a single one, to the 



energetically frozen out z degree of freedom while providing more states for xy degrees of 
freedom where correlated dynamics might be expected. The two-dimensional xy SPFs can 



then either be represented by primitive mode combination (figure 4(d) ) or by a multi-layer 



expansion, which can be favourable for bosons in traps being asymmetric with respect to 



the xy plain (figure 4(c)). Furthermore, the particle multi-layer expansion easily allows for 



simulating bosons with internal degrees of freedom. Figure 4(b) illustrates this with an 
example of one-dimensional, unpolarized spin-one bosons. As the z component of the spin 
is not assumed to be fixed, the spin degree of freedom S may occupy as many states as given 
time-dependent basis functions, namely three. 

After having introduced the species- and particle multi-layer MCTDHB wave function 
expansions, we briefly comment on two further tree diagram types: Naturally, the question 
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(c) (d) 

Figure 4: Examples for particle multi-layer MCTDHB: Figures a), c) and d) represent a 
species of bosons living in three-dimensional coordinate space. The wave function 
expansion a) reflects the equivalence of all three spatial dimensions, while in figures c) and 
d) the z-component is singled out by the fact that only one z orbital is provided. In figure 
c) , the xy SPFs are expanded according to the ML-MCTDH philosophy treating both the 
x and the y degree of freedom equivalently. In figure d), mode combination is applied for 

expanding the xy SPFs. A particle multi-layer expansion is used for describing 
one-dimensional spin-one bosons in figure b). The spin degree of freedom is denoted by S. 
The dashed lines in all figures shall indicate that the shown species branch could be a part 

of a species multi-layer expansion. 

arises whether tree diagrams are possible in which a bosonic node has bosonic children 
nodes for describing bosons made out of bosons. In the outlook, an extension of ML- 
MCTDHB to fermionic species will be discussed. Especially in this context, the question 
above gains importance as such trees could possibly mimic the forming of cooper pairs out 
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of electrons. Unfortunately, such tree structures break the permutation (anti-) symmetry 
of the constituting bosons (or fermions). Hence, this class of wave function expansion can 
only be applied as an approximation being valid under certain circumstances and will not 
be considered in this work. 

Finally, we would like to point out that the second layer in species ML-MCTDHB, i.e. the 
species layer, does not need to consist only of bosonic nodes alone but can also contain normal 
nodes representing additional distinguishable degrees of freedom. In such a setting, we would 
say that the bosonic species nodes span bosonic subbranches in the tree presentation, while 
the normal nodes of the second layer induce normal subbranches. In the following section, 
it will become clear how one can easily derive the equations of motion for the coefficients 
belonging to nodes of a normal and a bosonic subbranch. 

Summarizing, ML-MCTDHB makes use of the multi-layer concept above the bosonic 
nodes or below the child node of a bosonic node. Whereas the former can be used for an 
tailored treatment of inter-species correlation, the latter can be employed for reflecting the 
geometry of the system in coordinate space, i.e. the trap, and/or for considering internal 
degrees of freedom. 

2. Correspondence principle between ML-MCTDH and ML-MCTDHB 

The ML-MCTDHB equations of motion, which lead to a variationally optimal representa- 
tion of the total wave function at each instant, can be derived by applying the Dirac-Frenkel 
variation principle to the ML-MCTDHB wave function ansatz. From a technical point of 
view, however, it is much easier to study the relationship between ML-MCTDHB and ML- 
MCTDH firstly and to deduce the ML-MCTDHB equations of motion from ML-MCTDH 
thereafter: 

When dealing with systems containing indistinguishable bosons, a natural way of mode 
combination is to group all the bosons of the same species into one species node, i.e. one 
bosonic node, and this is the common starting point of ML-MCTDH and ML-MCTDHB. 
The difference between these two methods is that the SPFs of the species node are expanded 
with respect to permanents in ML-MCTDHB rather than Hartree products (see (T2B|) ). and 
consequently all identical bosons are represented by a single subnode below the species node. 
The correspondence principle now states the following: (i) Any ML-MCTDHB expansion can 
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be unravelled into a corresponding ML-MCTDH expansion, (ii) The expansion coefficients 
of non-bosonic nodes in ML-MCTDHB are the same as their counterparts in ML-MCTDH, 
and there is a linear mapping between the coefficients of a bosonic node in ML-MCTDHB 
and their ML-MCTDH counterparts. Consequently, we may translate the ingredients of 
the equations of motion of section III Al and the equations itself into the ML-MCTDHB 
framework. 

To demonstrate the mapping between the bosonic coefficients and their counterparts, 
we consider a subbranch of a ML-MCTDH expansion, which contains N identical bosons. 
Paying respect to their indistinguishability, we provide m orbitals for each boson and group 
them on the same layer in the ML-MCTDH tree. Such a combination of all the bosonic 
DOF to a single logical DOF, the species DOF, appears quite natural as the permutation 
symmetry of the bosons enforces strong correlations between them. 

Each species SPEfl, say \ijj r ) for simplifying the notation, is expanded in terms of Hartree 
products made of the iV SPFs. Each factor of such a Hartree product is an element of 
the set {\<Pi) \ i € [l,m]}. The indistinguishability of the bosons is then reflected in the 
symmetry of the expansion coefficients A r -i with respect to permutation within the multi- 
index / = (ii, ijf). As a consequence, the SPF \ip r ) can be equivalently expanded in 
permanents made of the orbitals \ipi): 

|Vv> = A r-j\fh) ■ ■ ■ \fi N ) = A-,fi\n)- (27) 

I n\N 

This identity implies that for each such a ML-MCTDH expansion of a bosonic species there is 
an equivalent ML-MCTDHB expansion and vice versa. An example for this correspondence 
principle is given in figure [5j Since we require the permanents to be orthonormal, i.e. 
(m|n) = Srn,rt, we arrive at the following correspondence relation between the ML-MCTDH 
and ML-MCTDHB coefficients A r -i and A r .ft, respectively: 



A r -n = \\yfm r A r -i, (28) 
1L=i n i- 

where / is some multi-index, which contains n s times the number s, s G [l,J7i]. This 
translation rule ensures the orthonormality of the species SPFs in both representations, i.e. 
'Ylii A* k .jA r -j = Skr if and only if A* k . n A r] fi = Sk r - Because of this one-to-one correspon- 
dence between the coefficients, we omit the tilde in the notation. The translation of the 



Or the total wave function in case of only a single species. 
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(a) (b) 



Figure 5: Visualization of the correspondence principle between ML-MCTDHB and 
ML-MCTDH: In subfigure a), a bosonic branch of a ML-MCTDHB wave function 
expansion is shown. This branch represents three bosons, which may occupy two orbitals. 

The dotted lines starting from the bosonic node as well as from its child node shall 
indicate that the shown branch may be incorporated in a species- or particle multi-layer 
MCTDHB expansion or both. Subfigure b) shows the corresponding ML-MCTDH 
representation of this branch. Please note that each of the three children nodes is 
associated with two SPFs because of the indistinguishability of these DOF. Due to the 
same argument, the very same subtree must be beneath each of these children nodes, 
which is again indicated by the two dotted lines. 

equations of motion from ML-MCTDH to ML-MCTDHB mainly subjects to the replace- 
ment of the summation over Hartree products of the species node to that over permanents. 
The replacement can be illustrated as: 

E /(-"«) = E fT^i/(\f^fr-^A < 29 > 

/ n Ai— 1 * 

where / is an arbitrary function. Following this, we can see that the equations of motion 
for the DOF above and below the bosonic node(s) are identical between ML-MCTDH and 
ML-MCTDHB. So the major difference between the two methods shows up in the bosonic 
node(s), where the density and the mean field (operator) matrices have to be expressed in the 
second quantization language as well as their action on the SPFs. In the following sections, 
we will apply the correspondence principle in order to derive the equations of motion for 
species- and particle ML-MCTDHB. 

Before we step to the concrete equations of motion of ML-MCTDHB, we would like to 
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comment briefly on conservation laws: As shown in the appendix of [8], all sets of equations 
of motion which can be derived from the Dirac-Frenkel variational principle preserve both 
the norm and the energy of the total wave function. In addition to this, we have proven 
analytically that the equations of motion of the MCTDH family preserve a certain class of 
symmetries in the following sense: Let the Hamiltonian H have some symmetry group Q 
and let G be an element of the unitary representation of Q in the many-body Hilbert space, 
i.e. [G, H] = 0. Here, we only consider symmetry operations which can be decomposed 
as G = ^ i c/i, where gi is an element of the unitary representation of Q in the Hilbert 
space of the z-th physical DOFj For simplicity, we focus on the MCTDH method^. If 
at t = the many-body wave function |^(0)) and all the SPFs |0j(O)) are now invariant 
under G and g~j, respectively, they remain invariant foreveio, i.e. G\^f(t)) = e* e |\l/(t)) and 
9j\4>j{t)) = el6j \ ( t ) j{^)) f° r a ^ t- This property of the MCTDH family cannot be expected a 
priori due to the truncation of the many-body Hilbert space and the complicated integro- 
differential equations for the SPFs and we are not aware of any work which has shown this 
before. Possibly this symmetry conservation can be employed for specifying the configuration 
selection in selected configuration MCTDH type methods (cf. 38| and references therein). 
For a sketch of the proof, see appendix |A] 



3. Species Multi-Layer MCTDHB 

In this section we present the theory of multi-layer MCTDHB in detail, and firstly we 
focus on species multi-layer, considering a mixture system of S species confined in one- 
dimension, as shown in figure 2. As indicated in equation ( |25|) . the total wave function is 
expanded as: 

Mi m s s 

i*> = £-£<....* ni«C>= £ ^i j1 >> (3°) 

h=l ig=l «=1 I=(ii,i2,...,is) 

where all time-dependencies are omitted as always. The species SPF \<pf K ) is in turn repre- 
sented in terms of permanents \n) K defined as in equation (Il7p : 

I0? K > = £4 ! »- ( 31 ) 

h\n k 

6 Some of the gj-th may be the unit operator. 

7 The sketched proof also holds for MCTDHB, MCTDHF and ML-MCTDH and thus also for ML-MCTDHB. 

8 In ML-MCTDH, the invariance of the SPFs on all layers under the corresponding symmetry operations 

has to be assumed at t = and will then be preserved forever. 
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The index k of the permanent shall indicate that \n) K is a permanent made of the SPFs 
which are associated with the physical DOF of species k bosons, i.e. x K . In the end, 
these SPFs are expanded with respect to the time-independent PBFs \x*): 

l^) = E4riO- (32) 

r 

Without loss of generality, we restrict ourselves to Hamiltonians involving at most two-body 

terms - a generalization to R > 2 body interactions is straightforward. So we assume a 
Hamiltonian of the form: 

s 

H = Y^\Hk + V k ] + J2 ( 33 ) 

where H K contains all one-body terms of species k, i.e. H K = Ei=i ^ with hf acting only 
on the z-th K boson. V K contains all intra-species interaction terms of species k and equals 
Yli<j Vij with vfj affecting only the i-th and j-th k boson. Finally, W KtK > denotes the inter- 
species interaction between the species k and k' and equals E^Ti ^2f=i with w*'* acting 
only on the i-th k and the j-th k' boson. Without loss of generality, all two-body operators 
are assumed to be given as a sum over products of two one-body operators: 

*s= E<^r«r> (34) 

V 

*■ k,k' \ v t^k.k' ~ KM a k'v 
W i,3 = l^ D v W i W j ■ 

V 

Such a form is particularly handy for avoiding high-dimensional integrals when calculating 
matrix elements and for the sake of a recursive formulation of the ML-MCTDHB ingredients. 

For applying the correspondence principle, we firstly summarize that the total wave 
function and each SPF is expressed as a summation over coefficients times states |C), which 
are either Hartree products (for the top layer), permanents (for the species layer) or PBFs 
(for the physical layer). According to the correspondence principle, the top and the non-top 
node equations of motion, respectively, then read: 

id t A l c = Y,(C\H\D)A l D) (35) 

D 

= E^K 1 " pa '^ E^T 1 * (HY'\, m Al« D \D), (36) 

D m 

where the density matrix, the mean field operator matrix and the projector are given by 

tively. Here, the SHFs are defined in the same way as in (Q everywhere in the tree but 
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for bosonic particle nodes, where we have the definition.]: 1^^) := -^==| l I / ). Consequently, 
we are left to evaluate those ingredients as well as the top layer mean field (C\H\D). In 
appendix [HI we give in detail the equations of motion for the coefficients of different layers 
as well as the formulas for all the relevant ingredients. 

4. Particle Multi-Layer MCTDHB 

In this section we turn to particle ML-MCTDHB, which is mainly designed for bosonic 
systems with a high-dimensional configuration space, spanned by the spatial DOF, and/or in- 
ternal DOF. To demonstrate the implementation of the particle-ML MCTDHB, we consider 
a single-species bosonic system of N identical bosons living in a Q-dimension configuration 
space, spanned by the physical DOF {xi,x%, ...,Xq}. In ML-MCTDHB, we firstly perform 
a mode combination of these physical DOF, and generate the logical DOF one layer above 
the physical layer, as in equation (jSj), and repeat the mode combination recursively layer 
by layer, until the layer containing only one node, which refers to the DOF of the boson 
x = (xi,X2, ...,Xq). Above this particle node, there is the top node representing the state 
vector of the whole system. This mode-combination procedure corresponds to the tree 
structure shown in figure [BJ where the dashed box between the particle node and the physical 
nodes indicates the variety of possible mode-combination schemes. As an example, we have 
shown different mode-combination schemes of two or three physical DOF into one particle 
DOF in figure HI 

We associate a set of SPFs to each logical DOF on the non-top layer, and expand them 
with respect to the SPFs of the corresponding DOF one layer below, as: 

m 

ht)= E ^ li2 ,..^nic +1;Kr )^E^i^ ;K }> (sr) 

11,12, —,im r=l I 

where is the n-th SPF of the K-th DOF on the a-th layer. The state vector \^f) of the 

top node is, however, expanded with respect to permanents : 

\*) = J2A 1 n\fi), (38) 

Fi 

9 Please note that although the SHF for a bosonic node is given in the second quantization representation, 
we exclusively treat problems with fixed particle numbers in the ML-MCTDHB formalism. This in turn 
means that projections like ( V E'^' K |\E') have to be understood as an incomplete scalar product in the N a 
particle Hilbert space resulting in a one-particle state and not as the scalar product in the Fock space. 
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Figure 6: Illustration of a general particle multi-layer MCTDHB wave function expansion 
for spin-one bosons in three-dimensional coordinate space. The dashed box indicates the 
various possibilities for mode combination of the physical DOF, i.e. the spatial DOF x, y, z 

and the spin projection on the z-axis. 



where \n) is the number state as defined in equation ( fTTI) . The physical DOF are also 
associated with a set of SPFs, and these SPFs are expanded with respect to time-independent 
PBFs as in ML-MCTDH (see equation ([5])). If the physical DOF $ K IS cL scalar quantity 
like the position in one direction, the associated SPFs are represented by one-dimensional 
PBFs. In the case of strong correlations treated with primitive mode combination, i.e. if 
the physical DOF x K is vector- valued, its SPFs are still expanded with respect to PBFs, but 
each PBF is then a product of one-dimensional time-independent basis vectors: 



4rn 



K,<J\ 



(39) 



i=(ii 



An example for incorporating primitive mode combination in particle multi-layer MCTDHB 



is shown in figure (4(d)), where the x and y DOF share the same set of SPFs, and are 
grouped into one primitive node. 



For simplicity, we restrict ourselves to two-body interactions and, moreover, assume the 
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product representation of the Hamiltonian as usually: 

H = Y,Hi) + Y,V{i,j) (40) 

i i<j 

i i<j v 

h(i) and V(i,j) are the one-body Hamiltonian acting on the i-th boson, and the two-particle 
interaction operator acting on the i-th and j-th bosons respectively, where the two-particle 
interaction operator is written as the summation of products of single-particle operators 
acting on the i-th and j-th boson separately. The single-particle operators h, v % v (i=l,2) 
can be further expressed as the summation of products of operators acting on each of the 
physical DOF, e.g.: 

where h K „ is a operator acting on the physical DOF x K . 

According to the correspondence principle, the equations of motion for the non-top nodes 
are equivalent to those in ML-MCTDH, and the top-node equations of motion are expressed 
in the second quantization picture. The equations of motion of all the coefficients can be 
given as following: 

id t Al = ^2(n\H\m)A^ (42) 

m 

idt^z = E^K 1 - /'- :K )/'r •/• i:k m;;::; 

+ E^{E( r;K i( 1 - J?a;K )^i ja;K )} 

V J 

For the equations of motion of the non-top node, P a ' K = |</>"' K ) (4>%' K \ is the projec- 
tion operator, and h a ' ,K and v"l* (i=l,2) are the single-particle operators of the one-body 
Hamiltonian and two-particle interactions, which are defined recursively as shown in the 
ML-MCTDH section (cf. equation ([8])). p a,K and (v^) denote the density and mean field 
matrices of the corresponding node, respectively. For the particle node, the density matrix 
is calculated as: 



P*3 



2;1 = Jf E ^ + i)K- + i)(4 + ^4 + i> (43) 



n\N-l 
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and the mean field matrices for the two-particle interaction terms are calculated as: 

(v 2 Ji)ij = 4 Yl S ^ + i)K + + 1) \/(^ + X )K + hi + 1) (44) 



n|7V-2 p,q 



*(^)*4 + ^(« 1 i^l« 1 )- 



The density and mean field matrices of the nodes under the particle node are generated 
as in ML-MCTDH (see equations (j!3|14p ). which is again related to the correspondence 
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principle and enables one to apply the well-known recursion schemes 

By now we have introduced the particle-ML MCTDHB with an example of a single-species 
bosonic system. The generalization to bosonic multispecies systems is straightforward, and 
we just need to apply the correspondence principle to the layers below the particle nodes, as 
we have done here for the single-species case, while we apply the species multi-layer to the 
layers above. This way we also complete the merging of the species multi-layer MCTDHB 
with particle multi-layer MCTDHB, obtaining the ability to deal with bosonic systems of 
arbitrary species in arbitrary high dimensions. 

As a final remark, we would like to point out the relationship between ML-MCTDHB and 
the standard MCTDHB. For single-species system, ML-MCTDHB and MCTDHB coincide 
for one- dimensional bosons without internal degrees of freedom. When considering bosons 
living in a higher-dimensional configuration space, ML-MCTDHB reduces to MCTDHB if 
primitive mode combination is applied directly below the particle node. For mixtures, the 
state vector of the total system is expanded in standard MCTDHB as the summation of 
the products of the permanents states belonging to each species, while in species multi-layer 
MCTDHB the state vector is firstly expanded with respect to Hartree products of species 
SPFs, and the species SPFs are then expanded within the permanents basis of each species. 
Hence, ML-MCTDHB recovers the standard MCTDHB treatment if as many species SPFs 
are provided as there are configurations for each species. 

III. APPLICATIONS 

Having introduced the theory of ML-MCTDHB, we now turn to some applications of 
;he method. Our ML-MCTDHB code is based on the general ML-MCTDH implementation 
in its recursive formulation 15 1. Since this ML-MCTDH implementation is capable of 



handing arbitrary ML-MCTDH wave function expansions, we have extended the code to 
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species and particle multi-layer MCTDHB as well as mixtures of both. Hence, in principle 
an arbitrary number of bosonic species in arbitrary dimensions can be simulated as long as 
the necessary number of configurations is not too large. 

At the moment, our ML-MCTDHB code can deal with arbitrary two-body interactions 
- an extension to three- or four-body interactions would be straightforward. As explained 
above, the recursive ML-MCTDH formulation demands a certain product representation 
of the interaction potential and so does our ML-MCTDHB implementation. This product 



representation can be achieved by the POTFIT algorithm |39l. |40| . Since most of our further 
applications are situated in the field of ultra-cold bosonic atoms, the contact interaction will 
be extremely relevant. Approximating the delta interaction potential by a narrow Gaussian, 



the POTFIT machinery can be applied for obtaining the desired product form (cf. 12]). 
In two- and three-dimensional systems, a word of caution is in order here: For instance, 
when expressing the regularized contact interaction in three dimensions, V(f = X\ — x<i) = 
g 5^ (r) d r r, in terms of the particle coordinates xi, x%, one also has to bring additional terms 
in the product form due to the modulus of r and its derivative with respect to the different 
components of x\ and x*i- For one-dimensional settings, we have also developed an exact 
implementation of the contact interaction V{x\ — x-i) = g8(x\ — x 2 ) without reshaping it into 
the product form in order to liberate the simulations from the artificial length scale induced 
by the Gaussian approximation of the delta function (cf. appendix [C]). The main difference 
in the derivation of the equations of motion lies in the fact that for a Hamiltonian being 
not of the product form (J6]) one cannot separate the mean field operators into a scalar mean 
field matrix and a operator valued factor anymore, which breaks the recursive multi-layer 
formulation. 

Depending on the concrete problem, we employ an appropriate discrete variable repre- 



sentation (DVR) of the SPFs and the Hamiltonian (45] . So in the case of the exact imple- 
mentation of the delta interaction, the grid point spacing defines the smallest length scale 
and, thus, an ultraviolet cutoff. 

For integrating the equations of motion, we employ either ZVODE, a variable-coefficient 
ordinary differential equation solver with fixed-leading-coeflicient implementation [46], or 
DOPRI, an explicit Runge-Kutta method of order five with stepsize control and dense output 
47|. 



In the following, we focus on the tunneling dynamics in a one-dimensional double well 
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system made of a harmonic trap superimposed with a Gaussian at the trap centre, i.e. 



Vtrap(%) — x 2 /2 + h/V^ns 2 exp(— x 2 /2s 2 ) in harmonic oscillator units h = m = oj = 
1. Firstly, we will show the correctness of our implementation by comparing the single 
component tunneling results with results from the well tested implementation of MCTDH 
J48 1 . Then we apply species ML-MCTDHB to simulate the tunneling of a binary mixture 



16] . Afterwards, we turn to a 



and compare with the established ML-MCTDH approach 
more complex scenario of three species with higher number of particles studying signatures 
of correlations and discussing the scaling in comparison to MCTDH, ML-MCTDH and 
MCTDHB. Applications of particle ML-MCTDHB will be presented elsewhere. 

In all examples, we prepare the initial state by modifying the trapping potential such that 
it becomes energetically favourable for the bosons to be in the left well. All bosons are then 
put into the ground state of the resulting single particle Hamiltonian and, afterwards, we 
let the many-body system relax to its ground state by propagating the equations of motion 
in imaginary time. The resulting many-body state is finally propagated in real time in the 
original double well trap. 



A. Single-Species Tunneling 

As a first application, we simulate the tunneling dynamics of a single-species bosonic 
ensemble in the double well. The tunneling dynamics of a bosonic ensemble in a one- 
dimensional double well has been extensively studied, for both microscopic and macroscopic 
systems, with boson numbers ranging from two to the order of 10 6 and even mor e |12 



52]. Experiments on the double well tunneling have also been carried out [53 
theoretical approaches have been employed, for instance, the Gross-Pitaevskii equation 



50] . the Bose-Hubbard model 



51 




5211 as well as ab-initio methods: MCTDHB has been 



applied to the many-body system [31], and calculations via MCTDH have been carried out 
for few-body systems jl^j]. The simulations based on the single-band approximation predict 
that the tunneling is suppressed for interaction strengths above some critical value, while 
the extended model predicts the weakening of such a suppression in the strong interaction 
regime, where higher bands effects cannot be neglected. In this way the double well potential 
manifests itself as a proper test bed for the higher band effects. In this section we will 
simulate the double well tunneling with ML-MCTDHB, and perform a detailed analysis to 
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Figure 7: The population oscillation of four bosons in the left and right well, with (a) 
g=0.0, (b) g=0.5, (c) g=2.0, and (d) g=4.0. For interaction strengths g < 2.0 (figures 
(a)-(c)), six orbitals are enough to achieve a good convergence, while in the case of g = 4.0 
(figure (d)), ten orbitals are needed to meet the convergence requirement. The results of 
ML-MCTDHB are shown as solid lines, and the MCTDH results as crosses. A good 
agreement is observed between the results of both methods. 



resolve the higher bands effects. 

Here we present the tunneling dynamics of four and ten bosons in a double well potential. 
A sin-DVR is employed, which intrinsically introduces hard-wall boundaries at both ends of 
the potential (cf. appendix of [8|). Firstly we show the simulation of the four-boson tun- 
neling, in comparison with the MCTDH simulation using the Heidelberg MCTDH package 



48|. To have a direct comparison between ML-MCTDHB and MCTDH results, we adopt 
in both simulations a narrow Gaussian interaction to model the contact interaction, and the 
interaction is written as V(x\ — £2) — di^^ 2 ) 1 ^ 2 exp(— (x\ — X2) 2 /(2a 2 )), with a = 0.05. 
Next we extend the simulation to the ten-boson case, and the contact interaction is mod- 
eled by the exact delta function. The simulation of the tunneling of ten bosons becomes 
impractical with MCTDH, and we only perform the simulation with ML-MCTDHB. 

Figure [7] summarizes the population evolution of 4 bosons in the double well at different 
interaction strengths. When the interaction strength is zero, the system undergoes Rabi 
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oscillations, as shown in figure [7(a). As the interaction strength increases to 0.5, the oscil- 
lations almost vanish on a relative long time scale, which indicates the delayed tunneling 
behaviour. When the interaction strength increases to 2.0, as shown in figure [7(c), the am- 
plitude of the population oscillations is increased from less than 0.06 in delayed tunneling 
(figure [7(b)) to around 1.0, which is referred to as enhanced tunneling. As the interaction 
increases even further to 4.0, the equilibration state emerges during tunneling, where the 
populations of the left and the right well approach the value of two with only small fluctu- 
ations, as shown in figure [7(d). Summarizing, figure [7J illustrates the tunneling transition 
from Rabi oscillations through delayed tunneling and enhanced tunneling to the equilibration 
state as the interaction strength increases from zero to the strong interaction regime. The 
ML-MCTDHB results show a very good agreement with the MCTDH calculations. Only for 
the interaction strength g = 4.0, deviations occur, which can be explained by the different 
implementations of the POTFIT algorithm in MCTDH and ML-MCTDHB. Nevertheless, 
we still observe qualitatively the same behaviour of the emergence of the equilibration state 
in both simulations. Moreover, more orbitals are needed to achieve good convergence in 
the strong interaction case of g = 4.0, and we supply ten orbitals to this case, where good 
convergence can be deduced from the natural populations discussed in the following. 

The ML-MCTDHB simulation supplies the dynamical information of the system in terms 
of the time propagation of the wave function, and we can extract various dynamical proper- 
ties of the system, via different analysis procedures. Among other things, we can calculate 
the eigenvalues and eigenvectors of the reduced one-body density matrix, the so-called natu- 
ral populations and natural orbitals, respectively, as well as the one-body density profiles at 
different times during the tunneling. The natural populations can confirm the convergence 
of the calculation and also manifest themselves as a measure of the fragmentation of the 
system. To uncover the physical effect giving rise to the enhanced tunneling and equilibra- 
tion state, we plot the natural populations and one-body density profiles of the four-boson 
ensemble at different time instants during the tunneling process in figures [7J Figure [8(a) and 
(c) show the natural populations and one-body densities for g = 2.0. In figure [5(a), firstly 
we see the lowest natural population saturates to a value less than 10~ 3 , which confirms the 
convergence of the simulation, and also that only two natural orbitals are macroscopically 
occupied. In figure [8(c), we observe that the profile in the left well remains as a Gaussian 
packet, while the profile in the right well presents a two-hump structure. The two- hump 
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Figure 8: The natural populations and one-body densities of the system of four bosons in 
the double well at different time instants, with (a) and (c) g=2.0, as well as (b) and (d) 

g=4.0. 

profile is a signature of the occupation of the first excited state in the right well, and this 
indicates that the enhanced tunneling is due to the higher band occupation, i.e. the in- 
terband tunneling. In the 2-boson ensemble [12] the enhanced tunneling only takes place 
in the fermionization regime, i.e. for the interaction strength approaching infinity, while as 
the number of bosons increases, it becomes easier to excite higher bands, and the enhanced 
tunneling arises even in an interaction regime far below the fermionization limit. 

The natural populations of g = 4.0, as shown in figure Efb) shows a good convergence 
of the simulation with the lowest natural population saturating well below 1 percent, and 
at g — 4.0 more natural orbitals contribute to the tunneling process, which suggests that 
fragmentation of the system and decoherence occur during tunneling. Figure [H^d) shows 
the one-body densities at different times for the interaction strength g = 4.0, where the 
equilibration state dominates the tunneling. During the tunneling, the one-body density 
profile presents multiple wiggles in both the left and right wells, which indicates multiple 
higher-band excitations in the tunneling process. In the exact quantum dynamics study of 
the double well system 31|, the equilibration state is explained by the decoherence of bosons 



between the left and right well, and the multi-excitation of higher energetic levels in the left 
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Figure 9: The population oscillation of ten bosons in the double well, with (a) g=0.5, (b) 
g=1.0. Enhanced tunneling is obtained in (a), where the amplitude of the population is 
around 2. The slow evolution to the equilibration state is obtained in (b). 

and right well suggests that the tunneling process involves a large number of higher band 
number states, and in consequence multiple tunneling channels. The dephasing between 
these tunneling channels can be a source of the decoherence of the bosons in the two wells. 

Figure [9] shows the tunneling evolution of the system of ten bosons in the double well 
with respect to the interaction strength. In the case of ten bosons, we focus on the enhanced 
tunneling and equilibration state for a sufficiently high interaction strength. The enhanced 
tunneling is observed with the interaction strength as weak as g = 0.5, as shown in figure 
[HI and at g = 1.0 the system slowly evolves to the equilibration state. Figures [TOT a) to (d) 
show the natural populations and one-body densities at different time instants during the 
tunneling process of figures |9^a) and (b), respectively. The convergence of the simulations is 
demonstrated by the lowest natural population, which saturates to the order of 0.1 percent 
and 1 percent in figures [TUl (a) and (b), respectively. In figure [T07 c). the two- hump profile 
in the right well indicates that the enhanced tunneling is again as a result of interband 
tunneling, and the multi-wiggle structure in figure [TOT d) suggests that multiple tunneling 
channels are involved and this can lead to the decoherence between the two wells, and 
consequently to the appearance of the equilibration state. 
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Figure 10: The natural population and one-body densities of ten bosons in the double well 
at different time, with (a) and (c) g=2.0, as well as (b) and (d) g=4.0. 



B. Mixture Tunneling 

Let us now study the tunneling dynamics of two and three bosonic species, which are 
loaded in the left well of the double well trap. All species shall have the same mass, which 
is set to one, and shall experience the same double well potential Vt rap (x). The intra-species 
interaction strengths g a of the contact interaction potentials V a (x) = g a 5(x), however, are 
assumed to be different for different species a. Furthermore, the inter-species interaction is 
also modelled by a pseudo-potential: V aa i(x a — x a i) = gaa'^ixa — x a >). All delta potentials 
are implemented numerically exactly as sketched above. We choose h = 3 and s = 0.2 for 
the height and width of the barrier, respectively, which leads to three bands each consisting 
of two single particle eigenstates. The lowest band is separated by an energy difference of 
1.63 from the first excited band and its level spacing amounts to AE w 0.23 leading to 
a tunneling period of T = 2tt/AE pa 27 for non-interacting particles. Assuming an equal 
number of particles of each species, N a = N, and not too different intra-species interaction 
strengths, we provide for each species the same number of species SPFs, M a = M, and 
particle SPFs, m a = m. For preparing the initial state of the mixture, we modify the double 
well trap by letting V trap (x) = 20.0 for x > and also close to the left end point of the 
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spatial grid. 

In order to illustrate the correctness of our implementation, we first consider a binary 
mixture made of 2 A and 2 B bosons. Here, we give an example where we have chosen = 
0.3, gs = 0.5 gA and gAB = 0.1 g&. Furthermore, we allow each particle to occupy m = 3 
particle SPFs and each species to be in M = 4 species states. In the following, we compare 
simulations of this system done with ML-MCTDHB and ML-MCTDH [16]. Although the 
initial state for the relaxation run is - from a mathematical perspective - perfectly symmetric 
with respect to particle exchange within each species, one has to pay attention to the initially 
unoccupied species SPFs in ML-MCTDH. These have to be symmetrized 'by hand' because 
otherwise the ML-MCTDH propagation will not preserve the exchange symmetry within 
each species. 

MCTDH and its derivative methods are proven to preserve both the norm and the total 
energy exactly (cf. appendix of Q). Using the ZVODE integrator with an absolute and 
relative tolerance of 10~ 10 and integrating 100 harmonic oscillator time units, the norm of 
the total wave function deviates from unity by a factor of 10~ 7 and the total energy is 
conserved up to a factor 10" 6 for both ML-MCTDH and ML-MCTDHB. 

Before discussing the convergence of the simulations with respect to the number of particle 
and species SPFs, we first summarize the results for different one- and two-body observables. 
Figure [TT1 shows the time evolution of the probability for finding an A (B) boson in the left 
well. One clearly sees that the tunneling period of the A bosons is enlarged in comparison to 
the Rabi-tunneling period, a manifestation of the well known delayed tunneling (cf. 12 



In contrast to this, the probability evolution of finding a B boson in the left well qualitatively 
resembles a first Rabi-cycle but afterwards also features delayed tunneling. This observation 
is quite plausible: Since both the gs and gAB are relatively small in comparison to gA, one 
expects that the B bosons require a longer interaction time in order to show an interaction 
induced effect. The impact of the different interaction strengths can also be seen in figure 
T2J While the probability for finding two A bosons in the same well is well above 0.5 
for most of the propagation time, the B bosons tend to stay in the same well less likely. 
This bunching tendency of the A bosons might be interpreted as the onset of correlated pair 
tunneling induced by the intra-species interaction (cf. |l2j). To the contrary, the probability 
for finding an A and a B boson in the same well fluctuates around 0.5 indicating that the 
bosons of each species tunnel independently. 
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Figure 11: 2 A and 2 B bosons initially loaded in the left well of a double well trap: The 
blue (red) line shows the probability for finding an A (B) boson in the left well versus 
time. Lines correspond to ML-MCTDHB results while crosses are obtained via 
ML-MCTDH calculations. Parameters are: = 0.3, gs = 0.5^ and g^B — 0.1 gA- Each 
species is allowed to occupy M = 4 and each boson m = 3 SPFs both in the ML-MCTDH 
and in the ML-MCTDHB simulation. The dashed vertical lines indicate the first three 
Rabi-tunneling periods of a single particle in this double well trap. 

For judging the convergence of the simulations, we plot the natural populations for dif- 
ferent subsystems. Figure [13] shows the natural populations corresponding to the reduced 
density matrix of the whole species A (or B). One clearly sees that after about 25 time units 
three species states contribute to the total wave function with weights of the order of 89%, 
10% and 1%. The fourth state contributes so little that it could be neglected without affect- 
ing the results. The natural populations of the reduced density matrix corresponding to an 
A and a B boson, respectively, are plotted in figure [HI Here again, we notice that the lowest 
natural population stays well below 1%, meaning that its natural orbital, the corresponding 
eigenstate of the respective reduced one-body density matrix, has only marginal influence 
on the result. Furthermore, we observe that two natural orbitals contribute with almost 
equal weight during certain time intervals. Hence, a mean field approximation would be 
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Figure 12: Evolution of several joint probabilities for the same setup as in figure dU The 
blue (red) solid line refers to the probability of finding both A (B) bosons in the same well. 
The probability of detecting an A and an B boson in the same well is plotted as the green 
solid line. All solid lines are obtained by ML-MCTDHB, while the crosses correspond to a 
ML-MCTDH simulation. Again, the first three Rabi-tunneling periods are represented by 

the dashed vertical lines. 



improper, which was to be expected for such a few-body system. Finally, we note that the 
ML-MCTDHB results excellently agree with the simulation performed with the well tested 
method ML-MCTDH. 

Let us now explore a more complex tunneling scenario of three species, A, B and C, 
each consisting of 6 bosons. This setup both unravels interesting correlation effects and 
illustrates the beneficial scaling of ML-MCTDHB by introducing the extra species layer. 
Here, we choose qa^Na — 1) = 0.2, gs = 0.75 qa and qab = 0.05 ^/QaQb- Furthermore, the C 
bosons are assumed to have no intra-species interaction, i.e. go = 0, but a weak attractive 
coupling to the bosons of species A and B: Qac — 9bc — ~ 0.5 g^B- 

From figure [T5J we see that the A, B and C bosons exhibit Rabi tunneling with respect 
to the tunneling period. The probability amplitude for the A and the B bosons, however, 
decreases in the course of time, which might be a sign of equilibration. This is also supported 
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Figure 13: The natural populations, i.e. eigenvalues of the reduced density matrix 
corresponding to the whole species A (or equivalently B) are plotted versus time for the 
same tunneling scenario as in figure [TTJ The solid lines refer to the M L-MCTDHB and the 
crosses to the corresponding ML-MCTDH calculation. 

by a low accuracy long time propagation up to 300 time units with m = M = 3, where we 
found that for t > 100 time units the probability amplitude for the A and the B bosons 
fluctuates around 0.5 with a fluctuation amplitude of less than 0.15 (data not shown). 
Figure [15] also shows the corresponding mean field results, which are obtained by setting 
m = M = 1. Then, the ML-MCTDHB equations of motion effectively reduce to a set of 
three coupled Gross-Pitaevskii equations. We clearly see that the decrease of the probability 
amplitude is a many-body property, not present in the mean field description. In contrast 
to this, the tunneling amplitude of the C bosons is not damped and its dynamics in the 
many-body calculation coincides with the mean field description. 

A further phenomenon not captured in the mean field picture is unravelled in figure 
[T6l The probabilities for finding two bosons of the same species in the same well oscillate 
between 0.5 and 1.0 with the frequency 2/T in the mean field calculation. In the many-body 
calculation, however, the probability for finding two A or B bosons in the same well features 
damped oscillations leading to a saturation of 0.73, which indicates a bunching tendency, 
while the probability of finding two C bosons stays oscillating between 0.5 and 1.0. 

In order to measure the correlations within a species, we compare the joint probability 
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Figure 14: The dynamics of the natural populations of the reduced density matrix 
corresponding to an A and a B boson are shown in (a) and (b), respectively. All 
parameters are chosen as in figure [TTJ The solid lines refer to the ML-MCTDHB and the 
crosses to the corresponding ML-MCTDH calculation. 



P(a, a) for finding two a bosons in the same well with the same joint probability in the case 
of statistical independence. If P(a; L) refers to the probability of finding a a boson in the 
left well, the latter joint probability reads Pi n d(a, a) = [P(er; L)} 2 + [1 — P(a; L)] 2 . Then the 
ratio f(a) := P(a,a)/P in d(cr,a) measures the intra-species correlation: f(a) equals one if 
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and only if the a bosons tunnel independently. From figure [171(a). we may infer that the 
ratio / stays one for the C bosons while for the A and B bosons / becomes larger than 
one and increases during the propagation time. Thus, the A as well as the B bosons show 
a bunching behaviour, which is not captured in the mean field picture, of course. Figure 
[TTT b) shows the ratio f(A, B) := P(A, B)/P inc i(A, B) corresponding to the joint probability 
of finding an A and a B boson in the same well. From this figure, we may infer that the 
presence of the weak attraction between the C bosons and the A (B) bosons reduces the 
bunching tendency of an A and a B boson being in the same well. 

For discussing the convergence of the simulation, the ML-M CTDHB calculation for m = 
M = 3 is compared with the results for m = M = 4 in figures [15] and [T6j Both the single 
particle and the two particle probabilities show an excellent agreement. Only for the joint 
probability of finding two particles of different kind in the same well, there are marginal 
deviations for long propagation times (not shown). Hence, we may regard the simulation 
as being converged. This judgement is also supported by the time evolution of the natural 



populations: From figure 18(a) , we infer that most of the time only two natural populations 
significantly contribute to the reduced density matrix of the whole species A (or B or C) 
and, hence, to the total wave function. For times larger than 69, a third species state gains 



a weight over 1%. Figure 18(b) shows that the initially fully condensed state of the A 
bosons evolves into a two-fold fragmented state. Increasing the number of particle SPFs 
from m = 3 to 4 just leads to a reshuffling of the third-highest natural population without 
affecting the results. The natural populations corresponding to an B boson show a similar 
behaviour due to the similar intra-species interaction strengths (not shown). In contrast to 
this, the C bosons stay in a condensed state. A (not fully converged) long time simulation 
with m = M = 3 shows that the state of the C bosons becomes depleted by 2% due to the 
weak inter-species attraction (also not shown). 

We finally compare the scaling of the methods MCTDH, ML-MCTDH, MCTDHB and 
ML-MCTDHB applied to this multi-species setup by investigating the memory consumption 
for storing the total wave function. This also provides us with a rough estimate for the 
performances of the methods. Let S be the number of species. If we denote the number of 
DVR grid points by n and neglect the symmetrization option as well as the possibility of 
primitive mode combination, m SN + SNmn coefficients have to be propagated in MCTDH. 
The ML-MCTDH expansion equivalent to our species ML-MCTDHB expansion requires 
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Figure 15: Tunneling dynamics of 6 A, 6 B and 6 C bosons, initially loaded into the left 
well, for the interaction strengths 9a{Na — 1) = 0.2, gs = 0.75 9c — 0-0, 
9ab = 0.05 y/g~A9B and gAc — 9bc — —0-^9ab'- The probabilities for finding an A boson, a 
B boson and a C boson in the left well are plotted as the blue, red and green line, 
respectively. The solid lines refer to a ML-MCTDHB calculation with three species and 
three particle SPFs per species, m = M = 3, while the corresponding crosses are obtained 

for m = M = 4. The dotted lines show the respective mean field calculation, i.e. 
m = M = 1. The first three Rabi-tunneling periods are represented by the dashed vertical 

lines. 

M s + S Mm N + S Nmn coefficients. For a direct MCTDHB expansion, one needs ( Af ^ m ] ^ 1 ) 5 + 
Smn coefficients. Finally, the species ML-MCTDHB ansatz consists of M s + SM( N ^~ l ) + 
Smn coefficients. Table |T] lists the memory consumption of the different methods for S = 3 
species, N = 6 bosons per species, n = 110 spatial gridpoints and varied numbers of m 
particle and M species SPFs. If the necessary number of species SPFs is not too large 
with respect to the total number of bosonic configurations, ( N ^ 1 ), ML-MCTDHB clearly 
requires less memory than all the other methods. In the case of the three species tunneling 
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Figure 16: For the setting of figure [T5| the time evolution of the probabilities for finding 
two A bosons, two B bosons and two C bosons in the same well are plotted as the blue, red 
and green line, respectively. Solid lines correspond to a ML-MCTDHB calculation with 
m = M = 3 and crosses to m = M = 4, while dotted lines are the respective mean field 

results. 



example and m = M = 4, the ML-MCTDHB wave function ansatz consists of two orders 
of magnitude less coefficients than the corresponding MCTDHB ansatz. For large numbers 
of species SPFs, however, a MCTDHB expansion becomes preferable with respect to the 
memory consumption, which can be seen from the m = 3 and M = ( N ^ 1 ) = 28 example 
in table DD 



IV. CONCLUSION AND OUTLOOK 

The methods MCTDHB and ML-MCTDH have different emphases: While MCTDHB 
aims at employing the bosonic particle exchange symmetry for obtaining a better perfor- 
mance, ML-MCTDH focuses on how to make use of system intrinsic correlations in order 
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Figure 17: For the setting of figure ITS'} the time evolution of /(er) (see definition in the 
text) is plotted for o = A, B, C as the blue, red and green line, respectively, in (a). Solid 

lines correspond to a ML-MCTDHB calculation with m = M = 3 and crosses to 
m = M = 4. For a mean field calculation, one finds f(a) = 1.0 (not plotted). The solid 
line in (b) shows the time evolution of the joint probability of finding an A and a B boson 
in the same well relative to that joint probability in the case of statistical independence. 
The dots correspond to the same quantity in the absence of species C. 

to reduce the complexity of the many-body wave function ansatz. In general however, the 
high flexibility of the ML-MCTDH wave function ansatz is incompatible with the generic 
correlations due to the bosonic exchange symmetry in a MCTDHB wave function expansion. 

In this work, we have derived and applied the novel ML-MCTDHB method, which over- 
comes the above conceptual problems and benefits from both the flexibility of a multi-layer 
structure of the wave function ansatz and from an explicit consideration of the bosonic ex- 
change symmetry. We have shown that such a symbiosis is indeed feasible for two - from a 
physical point of view quite natural - classes of wave function expansion schemes (and any 
combination thereof): (i) In species ML-MCTDHB, the total wave function of a bosonic 
multi-component system firstly is expanded in Hartree products made of states with each 
of them corresponding to a state of a whole species. Then each of these species states is 
expanded in permanents like in MCTDHB. (ii) In particle ML-MCTDHB, a single compo- 
nent bosonic system is considered with the bosons living in two or three dimensions and/or 
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Figure 18: The dynamics of the eigenvalues of the reduced density matrices corresponding 
to the whole species A (or B or C) and to a single A boson are shown in (a) and (b), 
respectively. The solid lines are the m = M = 3 ML-MCTDHB results, while crosses 
correspond to m = M = 4. The horizontal dashed line in (b) indicates an ordinate value of 
0.5, i.e. the perfect twofold fragmented state. All parameters are chosen as in figure T 



having internal degrees of freedom. There, the total wave function is expanded in terms of 
permanents and a ML-MCTDH expansion is applied to all the orbitals underlying those per- 
manents. Summarizing, the bosonic exchange symmetry is exactly and efficiently taken into 
account in ML-MCTDHB as any state of indistinguishable bosons is expanded in perma- 
nents. The multi-layer concept is then employed for obtaining an optimized wave function 
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Table I: The scaling of the methods MCTDH, ML-MCTDH, MCTDHB and 
ML-MCTDHB is compared for the case of three species, each consisting of 6 bosons. The 
number of particle SPFs, m, and (for the multi-layer methods) the number of species SPFs, 
M, are varied. The spatial grid is assumed to consist of n = 110 points. Each table entry 

contains the number of coefficients needed for the wave function expansion of the 
respective method and its ratio with respect to the number of ML-MCTDHB coefficients. 

ansatz guided by the correlations between different species (intra- versus inter-species corre- 
lations) or between different spatial directions (e.g. in quasi two- or quasi three-dimensional 
systems) or internal degrees of freedom. By means of a mixture tunneling scenario, we have 
illustrated the beneficial scaling of ML-MCTDHB. Like in all MCTDH type methods, the 
convergence of a simulation with respect to the number of provided states in the wave func- 
tion ansatz can be judged by the eigenvalue distributions of the reduced density matrices 
corresponding to different subsystems. 

The equations of motion of ML-MCTDHB can be derived by employing the Dirac-Frenkel 
variational principle ensuring that at any instant in time the propagated wave function 
remains variationally optimal within the class of wave functions given by the ansatz. Here, 
we presented an alternative but equivalent derivation, which uses the ML-MCTDH equations 
of motion as a starting point. Moreover, we have proven that all MCTDH type methods 
preserve an arbitrary symmetry of a given system in the sense explicated in section III C 21 
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as long as the corresponding symmetry operation can be written as a product of symmetry 
operations acting on only a single physical DOF. A generalization to genuine "many-body" 
symmetry operations acting on more than one physical DOF remains an open issue. 

We have implemented ML-MCTDHB on the basis of the ML-MCTDH code in such a 



general way that, in principle, we can deal with an arbitrary number of species in arbitrary 
dimensions - only limited by the number of states needed for a converged simulation, which 
depends on the system details, of course. Furthermore, the code can in principle also deal 
with hybrid systems such as a single ion coupled to an environment of indistinguishable 
bosons as long as all the interactions can be efficiently brought into the POTFIT product 
form. 

Furthermore, we have shown the first applications of the ML-MCTDHB method in terms 
of various correlated tunneling scenarios of a double well potential. Particularly for a setup 
with three species named A, B, C, where all interactions have been chosen to be repulsive 
except for the C bosons, which have no intra-species interaction but weakly attract bosons of 
the other two species, we have observed the emergence of intra- and inter-species correlations 
and signatures of equilibration: The probability of finding an boson of kind A (B) in the 
left well tends to equilibrate at a value of about 0.5, which is not captured in a mean field 
description at all. Similarly, the probability of finding two A (B) bosons in the same well 
tends to equilibrate at a value of 0.7. Such correlations can intrinsically not be described by 
the mean field ansatz. A more careful investigation has shown that the interactions lead to 
a bunching tendency of the A (B) bosons, which is increasing in time. We have shown that 
the weak attraction between the C and A (B) bosons reduces the slight bunching tendency 
of finding an A and a B boson in the same well. To summarize, there are various interaction 
induced correlation effects in such a three species tunneling setup, which can be studied 
systematically in further works. In particular, a discrimination of classical from quantum 
correlation appears to be a challenging perspective. 

Finally, we note that both conceptually and practically it is relatively straightforward to 
generalize the presented method to fermionic systems, i.e. to ML-MCTDHF, or to mixed 
bosonic-fermionic systems, i.e. to ML-MCTDHFB: According to the presented correspon- 
dence principle, one could start with an appropriate ML-MCTDH expansion in which all 
indistinguishable particles of one kind are grouped into a species node. The ML-MCTDH 
expansion of each of these species nodes would then have to be reformulated in the sec- 
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ond quantization picture, i.e. in terms of permanents or Slater determinants depending on 
whether the species is bosonic or fermionic. This approach is in particular compatible with 
the present ML-MCTDH code [16] and our ML-MCTDHB extension. 
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Appendix A: Proof of Symmetry Conservation in MCTDH and its Derivatives 

Here, we sketch the proof that MCTDH (and all its extensions) preserves both the sym- 
metry of the initial many-body wave function and the initial SPFs as explicated in section 
III C 21 According to a standard argument, the Hamiltonian matrix (I\H\J) in (ITT]) does 
not allow transitions from Hartree products |J) with G\J) = e* e '|J) to Hartree products 
1 1) with G\I) = e i@ "\I) and 6' mod 2vr ^ 6" mod 27r, when propagating the coefficients 
for a small time step At, i.e. (I\H\J) vanishes for such Hartree products. Hence, we are 
left to show that the SPFs preserve their symmetry up to 0(At 2 ) when propagating also 
them for a small time step. We note that the mean field operator matrices of a Hamil- 
tonian containing at most F-body terms can be expressed by the F-body reduced density 
matrices and one-body operators. If we group the initial SPFs \(j)j(0)) into classes of same 
6j mod 2tt (cf. section Hi C2j) . the F-body reduced density matrices and the inverses of the 
one-body reduced density matrices in (1T21 are initially block diagonal. With this, one can 
show igjdt\4>j(t)) = ie 10: >dt\4>j(t)) at t = 0. The same arguments hold for any later time steps 
as long as At is small enough and the symmetry conservation becomes exact for At — > 0. 
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Appendix B: Equations of motion for species-ML MCTDHB 



In this appendix, we demonstrate how to construct the equations of motion for the 
coefficients of each layer in species multi-layer MCTDHB. Firstly, let us state the recursive 
formulas for calculating the reduced density matrices in equations fl36|) . The density matrices 
of the species layer are computed in the well known MCTDH manner: 

£ = EK)^& (Bl) 

j k 

where J K indicates a summation over all species SPFs indices j a but j K . refers to the 
multi-index (ji, j K -i, n, j K +i, js)- For calculating the reduced density matrix on the 
physical layer, we have to evaluate the SHFs of a bosonic particle node at first: 

*?») := ^l*) = E l*f> ( E V(n, + 1)/N K |n) J , (B2) 

3 \n\N K -i J 



which leads to: 



M K 

k,l=l n\N K -l 



1 K * 

4' = F E # £ «« fe) A >'S + 3' < B3 > 



with Qn(i,j)= yj(ni + l)(n j + 1). 

In order to derive formulas for the mean field operators, we should also express the 
Hamiltonian (133]) in the second quantization language: 



"Ik 



4=E^% 3; "><A,;, (B4) 



2 



^ = E E E ^r^rit&^i^^i^i^a^a;^. 

^ p,g=l «,j=l 

Here, we have introduced the one-body operators /i K = ft,* an d ^ n,v = '^i' 1 ' to simplify the 
notation. Please note that the instantaneous matrix elements occurring in (]B4j) have to be 
evaluated by means of the PBFs according to (132]) . 

Going from the physical layer upwards, we construct the mean field operator matrices and 
state the equations of motion layer by layer. Due to the presence of the projector P a,K and 
the projection on |C) in (136]) . we have to calculate the matrix elements (F\ \E), 
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where \E) and \F) are two PBFs or permanents for the physical (a = 3) and the species 
(a = 2) layer, respectively. We are allowed to treat the one-body, intra-species and inter- 
species interaction terms separately. Furthermore, we may employ the factorization of the 
mean field operator matrix into an operator and a scalar prefactor, the mean field matrix, 
( TT4"|) due to the product structure (134ft . 

First, we investigate the physical layer. But before diving into the concrete calculations, 
a word of caution is in order here: In (IB2j) . we have introduced the SHF of a bosonic physical 
DOF x K , which clearly contains N K — 1 bosons of species k. However, the mean field operator 
(^^ K \H\^^) is in fact an operator acting on single particle states. So one must not plug 
in the Hamiltonian in second quantization representation (IB4[) directly, rather one should 
use the first quantization representation and reformulate the result in terms the second 
quantization language according to the correspondence principle (cf. also footnote [9]). 

The peculiarity of the one-body terms H K is that they turn out to be uncorrelated in the 
sense that (\l/^ ;K |-ff K |\l/^ K ) equals h K p^ plus scalar terms which are killed by the projector. 
Hence, the scalar factor of the surviving term will cancel the inverse density matrix in (IHoD . 
In contrast to this, the mean field operators of the interaction terms show non-trivial mean 
field matrices: 

W \'„ * 3 J?) = S VK (N K - I) ^^^I^W*) ^ + scalar terms, 

V 

(*' i iw^w«> = x;^r(^wi^i*^> + ( B5 ) 

V 

+ 5 aK N v (^ K \w v ' u \^ K }) w K ' u + scalar terms. 

We again do not have to care about the scalar terms since they will be killed by the projec- 
tor. Probably the easiest way for arriving at the above expressions is to think in the first 
quantization picture and to interpret the SPF \(j)f K ) of the bosonic physical DOF % k as a 
SPF of the "first" k boson in an - according to the correspondence principle - equivalent 
ML-MCTDH expansion of the total wave function. Then the SHFs \^ 3 k ' K ) describe states 
corresponding to the "second" up to the "iV K -th" boson of the species k and all the bosons 
of the other species. 

Now we may interpret the operators in the mean field matrices as one-body operators in 
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second quantization language, which enables us to apply the relation (jB2[) directly: 

= E^i^i^)(*n ; i4A^r), (B6) 

So we are left to evaluate (^ ;re |aj ci a K j| 1 I^ re ) and ( 1 I r ^' re |ap. i a< T) j| 1 i r ^ re ), which are essentially 
the reduced two-body intra-species density matrix of species k and the reduced two-body 
inter-species density matrix of species k and a, respectively. Thus, those quantities can also 
be employed for a physical interpretation of a ML-MCTDHB simulation. 



N K 

r ' s l\N K -2 



x A 2 ' K r . , P r (m,^ 



r,s,u,v l\N K -l 
X 



0|JVa-l 

where Pf{n,i) = a/ (l n + + iJPi + 1) an d denoting the reduced two-body inter- 

species density matrix in the respective species SPF representation: 

jna 

Here, J KfT refers to a summation over all indices but those of k and a and denotes the 
multi-index J Kcr with k being in the r-th and a being u— th species SPF. Hence, we arrive 
at the final result for the equations of motion of the physical layer: 

id t A% ^E^K 1 - P ^\<) A t+ (B9) 

+EEwi(i - p 3 -' K K' v K) E^ 3;K )- (^)^^S 

s n,m 

+ E E E - p 3iK K>io E(p 3;K )r™ (toss ^S- 

s (t^k v n,m 

Please note that the projector P 3;K = £V 10^X0^1 a l so has to be expressed in terms of the 
PBFs by means of (132|) . Furthermore, we have assumed that = D*>° without loss of 
generality. 
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Next we turn to the ingredients of the equations of motion for the species layer. In 
this layer, the SHFs \^ ,K ) describe states of all species but k. Following the same logic 
as above, the one-body and intra-species interaction terms, H K and V K , respectively, are 
uncorrelated terms on the species layer. Thus, we are left to determine the mean field 
operators corresponding to W^, which are only non-trivial for either rj or a equal to k. 
Here, we evaluate the mean field operators for rj = k after separating them into mean field 
matrices and operator valued terms (cf. §U 



(^IW^J?) = J2 D » a WOnm (BIO) 



By employing the second quantization representation (JB4J), we can state the operator valued 
factor: 

^ ;k = E^i^i^KA* ( B11 ) 

p,q=l 

as well as the mean field matrices: 

(^r>n ; m = E (^K'l^f) <*£1<^j W">- ( B12 ) 

Inserting the definition of the SHF leads to the following: 

«14Aj«"> = Efe E QdU) (•>;;;• y-C- r ( B13 ) 

r > s T\N a -l 

So the equations of motion for the species layer read: 

^4* = E W - p2]K ^ + ^)l™> A S+ ( B14 ) 

m|AT K 

+ EE^ ,CT E ("K 1 - **>?v> E^)- 1 (c^r^S- 

Again, the projector P 2;K has to be explicated by means of (131 p . Besides, one has to evaluate 
the action of the operators H K , V K and w^ K on the species SPF |0 2;K ). This is straightforward 
to do when employing the second quantization representations (IB4j) . (1B1 1|) and, here, we 



only state how the archetypes of a one- and of a two-bodjo operator act on such a species 



10 Here, we refer to only intra-species interaction terms as two-body operators in contrast to inter-species 
interaction terms, which we regard as products of intra-species one-body operators (cf. (|B4I) ). 
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SPF: 



alAM K ) = E GKm) l f+ (bis) 

When plugging (IB4|) and (lB15j) into (IB14[) . we arrive at the final equations of motion. Finally, 



we turn to the top layer equation of motion: 

M K 



i ^i=EE^iA + ^M' K ) 4+ (big) 

K p=l 

M K Ma- 

+EEE(Ci(Ci^i« K )i« CT ) 

k<ct p=l g=l 

The occurring mean field matrices can be calculated with the help of (IB15P and the Hamil- 



tonian in second quantization language (1B4|) . Here, we only state the result for the mean 



field matrix corresponding to the intra-species interaction: 

-, rn K 

efeia^ E E^^fi^i^^fi^i^f) (si?) 

L v 



2 



Appendix C: Comments on the Implementation 

In several of the equations above, there are summations over bosonic number state con- 



figurations. In a computer program, these sums can be 



performed convenient 



mathematical framework of combinatorial number systems 57| as pointed out in 



y in the 



25|. The 



idea is to map a number state n to a single integer, which serves as a unique address for 
picking the corresponding coefficient, say A?!£, out of an array. In 25], a one- or two-body 



operator is applied to a bosonic state as follows: A nested loop over all occupation numbers 
ni,...,n m is performed with the constraint that Yli n i equals the number of bosons (of one 
kind). Applying a one- or two-body operator results in removing and adding bosons in 
certain orbitals. The "scattered" number state can then be mapped to a single integer and, 
thus, the action of a one- or two-body operator reduces to re-addressing the number states 
(and multiplying the respective coefficient with a known prefactor). 
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In our ML-MCTDHB code, we employ a scheme which is different in two aspects: First 
of all, we would like to note that we have reduced the summation over all iV particle number 
states to a summation over all N — 1 (N — 2) particle configurations for all one- (two-) body 
operator related quantities. This reduction improves the performance due to the binomial 
scaling of the number of configurations with the number of bosons. Secondly, our code 
calculates several index mapping tables before the actual simulation: (i) There is a table 
from which, for every integer address / of a number state, the corresponding occupation 
number ni,...,n m can be read off. Hence we can avoid a nested looping over the occupation 
numbers rij by simply looping over the integer addresses I. (ii) For any integer address / 
corresponding to a N—l particle number state I, we store the integer address I' corresponding 
to the N particle number state I + i for any i — 1, ...,m in another table, (iii) Similarly, 
for any integer address I corresponding to an N — 2 particle number state I, we store the 
integer address I' corresponding to the N particle number state l + i + j with i,j = 1, .., m. 
As the formulas like ( 1B15[) indicate, applying a one- or two-body operator to a bosonic state 
effectively means finding the correct integer addresses of the number states from which and 
to which the annihilation and creation operators scatter. The tables (ii) and (iii) allow for 
a extremely fast and efficient identification of these addresses, while table (i) is used for 
reading off the underlying occupation numbers needed for the prefactors. 

Finally, we would like to comment about interaction potentials which are not given in 



the POTFIT form. Any potential can be brought into POTFIT form [39|, |40j. In general, 
such a representation requires as many POTFIT terms as there are gridpoints but often only 
a few terms are really needed for a good approximation of the potential, which improves 
the performance. If, however, many POTFIT terms are needed, this can slow down the 
simulation and it can be beneficial to work with evolution equations which are not based on 
a potential of the POTFIT form. For such a potential, the equations of motion can be de- 
rived by starting with (1B9|B14|B16P and then performing the summation over the POTFIT 



terms first and foremost, which recovers the interaction potential in its "original" form. As 
a consequence, one cannot split the mean field operators into scalar valued mean fields and 
operator valued factors anymore, which breaks the recursive structure of a multi-layer expan- 
sion and leads to more complicated objects. In particular, this procedure is straightforward 
for the contact potential in one dimension: After inserting the exact POTFIT represen- 
tation S(xi — x 2 ) = J dy5(xi — y)S(x2 — y) into the equations of motion, the summation 

52 



over the POTFIT terms, i.e. the integral over y, can be easily performed and one arrives 
at the ML-MCTDHB equations of motion without the product form assumption. We have 
observed a speed up by one order of magnitude and sometimes even more for this efficient 
implementation of the contact interaction in comparison to the equations of motion with a 
POTFIT representation containing as many terms as there are PBFs. 
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